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Abstract 

A new class of micromechanically motivated chain network models for soft biological tissues is presented. 
On the microlevel, it is based on the statistics of long chain molecules. A wormlike chain model is applied 
to capture the behavior of the collagen microfibrils. On the macrolevel, the network of collagen chains is 
represented by a transversely isotropic eight chain unit cell introducing one characteristic material axis. 
Biomechanically induced remodeling is captured by allowing for a continuous reorientation of the predomi- 
nant unit cell axis driven by a biomechanical stimulus. To this end, we adopt the gradual alignment of the 
unit cell axis with the direction of maximum principal strain. The evolution of the unit cell axis' orientation 
is governed by a first-order rate equation. For the temporal discretization of the remodeling rate equation, 
we suggest an exponential update scheme of Euler-Rodrigues type. For the spatial discretization, a finite 
element strategy is applied which introduces the current individual cell orientation as an internal variable on 
the integration point level. Selected model problems are analyzed to illustrate the basic features of the new 
model. Finally, the presented approach is applied to the biomechanically relevant boundary value problem 
of an in vitro engineered functional tendon construct. 
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1 Introduction 



Collagen is a fibrous protein secreted by connective tissue cells. The distinguishing feature of a typical colla- 
gen molecule is its long, stiff, triple-stranded helical structure, in which three collagen polypeptide chains, the 
so-called a chains, are wound around one another in a rope-like superhelix. So far, about 25 distinct collagen 
a chains have been identified. The a chains that make up type I collagen are by far the most common. Type 
I collagen is a fibril-forming collagen which is present in nearly all connective tissues such as bone, skin, ten- 
dons or ligaments. After being secreted into the extracellular space, collagen molecules assemble into collagen 
microfibrils. These are sparsely cross-linked higher order polymers which are about 10-300 nm in diameter and 
several hundreds of micrometers long, see e.g. the fundamental textbook on cell biology by Alberts et al. [1] 
or the rather mechanically oriented textbooks on soft biological tissues by Fung [18], Cowin & Humphrey [11], 
Holzapfel & Ogden [22], Humphrey [23] and Humphrey & Delangc [24]. 

Polymer chains have many conformations of nearly equal energy Perturbing the chains away from their equi- 
librium conformations typically generates entropic forces that oppose these perturbations. This is the basis for 
entropy based elasticity. Since the number of different configurations which a long chain molecule may assume 
is very large, the treatment of each of them individually would be a complex, maybe even unmanagable, task. 
Long chain molecules are thus commonly described by statistical mechanics, a concept which was originally de- 
veloped in the context of entropic rubber elasticity by Kuhn [28,29] or Kuhn & Grim [30], see also the textbooks 
by Flory [14] or Treloar [43]. As pointed out by Boyce [6] and Boyce & Arruda [7], statistical, microscopically 
motivated models are typically characterized by a limited number of well-defined, physically-motivated material 
parameters. This is in contrast to the macroscopic phenomenological models documented e.g. by Ogden [37,38], 
Treloar [43,44] or Holzapfel [21]. 

Collagen chains are extremely rich in proline and glycine. While the former stabilizes the helical conformation 
of each a chain, the latter allows the three a chains to pack tightly together in the final collagen superhelix. 
Unlike polymer chains in rubber, which are of rather uncorrelated nature, collagen chains in biological tissues 
thus have to be classified as correlated chains from a statistical point of view. Rubber chains are typically char- 
acterized as freely jointed chains, the configuration of which resembles a random walk, whereas biological chains 
correspond to "wormlike chains" with a smoothly varying curvature, see Kratky & Porod [26]. Traditionally, 
the wormlike chain model has been applied to describe the behavior of the DNA double helix, e.g. Marko & 
Siggia [34] or Bustamante et al. [8]. Only recently, the wormlike chain approach has been adopted to simulate 
the constitutive behavior of the collagen triple helix by Bischoff et al. [4,5] and Garikipati et al. [19]. 
From the standpoint of polymer structures, rubber as well as soft biological tissue consists of a complex three- 
dimensional network of polymer chains laterally attached to one another at occasional points along their lengths. 
To account not only for the behavior of the individual chains as such, but also for the characteristic cross-link 
effects of the network structure, a number of different chain network models have been proposed over the past 
60 years. The common feature of all these network models is a characteristic unit cell which is assumed to 
provide an adequate representation of the underlying macromolecular network structure. The existing chain 
network models can basically be classified in two categories: affine and non-affine models. 
The first affine network model was the three chain model by James & Guth [25] and Wang & Guth [47] . Relat- 
ing the overall strain state in an affine manner to the stretches of three mutually orthogonal chains, the model 
generally overestimates the overall stiffness when assuming that one out of these three chains is aligned with 
the direction of maximum principal strain. Moreover, due to the orthogonal arrangement of the three chains, 
cross-linking network effects arc basically neglected. To remedy these deficiencies, an affine full network model 
was suggested by Treloar & Riding [45], compare also Wu & van der Giessen [48] or Miehe [36]. As its chains 
are distributed with equal probability over the solid angle, the affine full network model is considerably weaker 
than the three chain model since only a limited number of chains are strained up to the locking stretch upon 
uniaxial deformation. 

The first representative of the class of non-affine models was the four chain tetrahedron model by Flory & 
Rchner [17] and Treloar [42]. It is not surprising though, that due to the tetrahedral shape of its unit cell, 
the four chain model reveals a non-physical anisotropic response. The non-affine eight chain model based on a 
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cubic unit cell is maybe the most prominent isotropic chain network model today, compare Arruda & Boyce [2] 
or Boyce [6]. Although counterintuitive at first glance, non-affine chain network models seem to be superior 
over afiine models in predicting the overall response under various load cases as they assume an instantaneous 
orientation of the unit cell with respect to the principal stretch space. Critical comparisons of the different 
isotropic chain network models can be found in the classical textbook by Treloar [43] or in the monographs by 
Flory [15], Boyce & Arruda [7] or Miehe et al. [36]. 

It was already pointed out by Kuhn & Grim [30] that, for rubberlike materials, the assumption of isotropy is 
typically only met in the small strain regime. Upon further loading, the polymer chains were found to reorient 
themselves with respect to the loading direction. There is experimental evidence that this effect of anisotropy 
is even more pronounced in soft biological tissues, see e.g. the early experiments by Lanir & Fung [32], or the 
textbooks by Fung [18], Holzapfel & Ogden [22] and Humphrey [23]. It is widely accepted that this effect of 
anisotropy can be attributed to the collageneous network structure in which the polymer chains are typically 
oriented with respect to the predominant loading direction. A classical example is provided by Langer's lines 
in skin, pointing in the direction of maximum tensile strain. Similar effects can obviously be observed in mus- 
cles, tendons and ligaments. A first attempt at modeling the anisotropic response of soft tissues based on an 
anisotropic eight chain unit cell has been presented recently by Bischoff ct al. [3-5] and Garikipati et al. [19]. 
In the present work, we shall adopt the above concepts to derive a transversely isotropic chain network model, 
see also Kuhl et al. [27]. Following the ideas of the classical eight chain model, the isotropic in-plane response 
is represented in a non-affine manner, while the out-of-plane stretch is related to the overall macroscopic strain 
through an affine transformation. 

Although characteristic for fully developed tissue, the pronounced orientation of collagen chains is typically 
not present in the embryonic state, see e.g. Calve et al. [9]. It is only upon mechanical loading, that the 
tissue exhibits a directional strengthening due the local rearrangement of the collagen fibers, an effect with 
is referred to as "functional adaptation" or "remodeling" in the biomechanical literature, see e.g. Taber [41]. 
Again, two approaches can be distinguished, macroscopic phenomenological concepts and micromechanically 
motivated strategies. The former are typically based on the introduction of a ficticious growth configuration and 
the characterization of its evolution with respect to the undeformed configuration, see e.g. Rodriguez et al. [39], 
Lubarda & Hoger [33], Epstein & Maugin [13] and Garikipati et al. [20]. Micromechanically motivated remod- 
eling theories, on the contrary, are based on the rigorous reorientation of collagen fibers, see e.g. Cowin [10], 
Vianello [46], Sgarra & Vianello [40], Driessen et al. [12]. The present work is particularly related to the recent 
contribution by Menzel [35] , who suggests a reorientation of the material's principal axis with respect to the 
predominant loading direction. 

Within this contribution, we suggest a gradual alignment of the axis of transverse isotropy with respect to the 
direction of maximum principal strain. Recall, that upon loading, the non-affine eight chain model by Arruda & 
Boyce [2] tacitly assumes an instantaneous reorientation of the representative unit cell, such that its axes always 
remain aligned with the principal strain axes. Herein, the approach of instantaneous alignment is adopted for 
the isotropic non-affine in-plane response, while for the affine out-of-plane response, we postulate a gradual 
reorientation. 

The present manuscript is organized as follows. Section introduces the basic ideas of the micromechanics 
of long chain molecules to simulate the individual collagen fibrils. In particular, two different types of model 
chains are discussed and compared: the freely jointed chain and the wormlike chain. Section [3] then deals with 
the mechanics of the chain network. Based on the concept of a representative eight chain unit cell, we consider 
a transversely isotropic network model of wormlike chain type. We then address the issue of reorientation. Ac- 
cordingly, section 01 focuses on the continuum model of remodeling and its algorithmic realization. The features 
of the single chain model, the chain network model and the reorientation model are illustrated individually in 
terms of a homogeneous model problem under uniaxial tension in each section. Finally, the overall model is 
used to predict the fiber reorientation of randomly oriented collagen fibers in a cylindrical model tendon. The 
suggested remodeling approach is discussed in section [5] 

Remark 1 (Notion of affinity) Throughout this paper, we shall apply the notion of "affinity" in the context 
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of "affine motion". According to Holzapfel [21], affinity in this sense implies that changes in the length and 
orientation of lines marked on chains in a network are identical to changes in lines marked on the corresponding 
dimensions of the macroscopic sample. In the principal strain space, the eight chain model can thus be classified 
as an affine model. For arbitrary load cases, however, the eight chain model does not assume affine deformation 
of all chains, see Boyce & Arrunda [7] or Miehe et al. [36]. Note that we do not use the notion of affinity in the 
context of "affine junction motion" which it was related to earlier in the context of constrained junction theories 
or phantom models, i.e. models which nicely caputre the network behavior at moderate strains, see e.g. Flory & 
Erman [16]. 

2 Micromechanics of a single collagen chain 

We shall begin by introducing the kinematics of a single polymer model chain. In the simplest case, this chain 
can be characterized through N rigid bonds of equal length I, the so-called Kuhn length, see e.g. Kuhn [28,29]. 
Accordingly, the total stretched-out length of a chain, the contour length, is given as L = N I. The deformation 
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Figure 1: Kinematics of individual chain models 



of this model chain is then typically defined in terms of the end-to-end distance r, i.e. the length of the vector 
pointing from one end of the chain to the other, whereby < r < L. Alternatively, the deformation of the 
chain can be characterized through the relative chain stretch A = r/L, i.e. the dimcnsionlcss ratio between the 
end-to-end length and the contour length. While I, N and L are constant for a particular chain, the end-to-end 
length r and the chain stretch A = r/L change as the chain is subjected to applied forces. Two different types 
of chain models can be distinguished from a statistical point of view: uncorrelated and correlated chains, as 
illustrated in figure ^ I n what follows, we will discuss different representatives of each of these classes, i.e. the 
classical single parameter freely jointed chain and the two parameter wormlike chain. 



2.1 Uncorrelated chains: The freely jointed chain model 

The most common model for chains is the random flight or rather freely jointed chain model. The freely jointed 
chain consists of N bonds of fixed bond length I, whereby the directions of neighboring bonds are completely 
uncorrelated in the sense that all directions for a given bond are of equal probability, irrespective of the directions 
of the neighboring segments. The model is thus characterized through one single parameter, the contour length 
L = N I. Figure^ lsffc, illustrates a freely jointed chain with TV = 20 bonds. Let us introduce the probability 
density p(A) , i.e. the probability that a chain of the contour length L takes a configuration characterized through 
the end-to-end length r. According to the classical Boltzmann equation s^ c — k ln(p), the entropy s^ c of a 
single chain can be expressed in terms of the Boltzmann constant k = 1.3810~ 23 J/K and the probability density 
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p. For purely entropic chains, the free energy tjj^ c of a single chain can thus be expressed as ip i]c = — k9\n (p) 
where 9 is the absolute temperature. Depending on the range of stretching, either Gaussian or non-Gaussian 
statistics are commonly applied in order to specify the particular entropy changes upon deformation. For the 
classical Gaussian case for which p — poexp(—3/2Nr 2 /L 2 ), the free energy of an individual chain takes the 
following form 



"r gau 



(1) 



where 1$ c is the value of the chain energy in the unperturbed state. The single chain force follows straightfor- 
wardly as 



f£ u = k6NZ- 



(2) 



and is thus a linear function of the relative stretch r/L. Alternatively, we could apply non-Gaussian 
statistics of inverse Langevin type as introduced by Kuhn & Grim [30]. With p = p exp(—N£~ 1 r/L — 
iVln^" 1 / sinh(£~ 1 ))), the strain energy ip^ c of the freely jointed chain can be expressed as 
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(3) 



where i/Jq C is the value of the chain energy in the unperturbed state and C 1 is the inverse Langevin function 
with C (r/L) = coth(r/L) — L/r. Accordingly, 



ft = k6NL- 1 



(4) 



defines the force stretch relation for the freely jointed chain with inverse Langevin statistics. Recall that the 



inverse Langevin function can be evaluated by a Pade approximation as C 



(3 - r 2 /L 2 )/(l - r 2 /L 2 )r/L, 



as in Miehe et al. [36]. At small stretches r/L for which C 1 « 3r/L the force of the inverse Langevin chain 
/lan = k6 NC~ l thus obviously approaches the force of the Gaussian chain /g J a c u = k6 N 3r/L . 



2.2 Correlated chains: The wormlike chain model 

The distinguishing feature of the wormlike chain or Kratky-Porod model chain is the continuity of the direction 
of its contour in space, see Kratky & Porod [26] or Flory [14]. As such, it is characterized through a smooth 
curvature whose direction changes randomly but in a continuous manner. This property is essentially reflected 
through a second parameter besides the contour length L = N I, namely the persistence length A. The persis- 
tence length can be understood as the sum of the average projection of all bonds onto the direction of the first 
bond, as sketched in figure Q right. Varying between / < A < L, the persistence length is thus a particular 
measure of stiffness, see also Landau & Lifshitz [31]. Accordingly, the persistence length of the uncorrelated 
freely jointed chain is equal to the length of the first bond A — I while the persistence length of an infinitely 
stiff chain with almost beam-like properties is equal to its contour length A = L. The strain energy ip wlc of the 
wormlike chain model 
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can be derived straightforwardly by integrating the force stretch relation for a wormlike chain 
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(5) 



(6) 



as originally suggested for the DNA double helix by Marko & Siggia [34] and Bustamante et al. [8] and applied 
for the collagen triple helix by Bischoff et al. [5]. Again, i/j™ 1c is the value of the chain energy in the unperturbed 
state. 
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2.3 Example: Comparison of different chain models 

To illustrate the fundamental differences between the freely jointed chain model and the wormlike chain model, 
we plot the different force elongation curves for two individual model chains. Figure [21 left, depicts the force 
elongation behavior of both, a Gaussian and an inverse Langevin freely jointed chain with the force / JC being 
scaled by the factor l/[k N] . The two curves clearly monitor the deviation of Gaussian and the inverse Langevin 
statistics in the large strain regime, for which the linear force elongation behavior of Gaussian statistics is no 
longer appropriate. The inverse Langevin freely jointed chain model, however, nicely captures the characteristic 
locking behavior close to the locking stretch r = L. In the low strain region, the chain shows nearly no resistance 
to loading while close to its full extension, the chain stiffness increases considerably. 



single chain force f - freely jointed chain model single chain force f - worm like chain model 




Figure 2: Force vs. elongation response of individual chain models 



Next, we elaborate the force elongation response of a single wormlike chain model. Figure right, shows the 
chain force / wlc scaled by l/[k 9] for varying chain stretches r/L at different persistence lengths, A, varying from 
A = 0.125, right curve, to A = 1.000, left curve. Again, the characteristic locking behavior is nicely captured by 
the model as r/L approaches one. However, in contrast to the single parameter freely jointed chain model, the 
wormlike chain model offers the additional freedom of a second parameter, namely the persistence length. With 
this second parameter, the wormlike chain model is not only able fit the locking stretch but also to capture 
the shape of the force elongation curve appropriately. For larger values of of the persistence length indicating 
initially stiffer chains, the locking response, i.e. the behavior in the r / L — ► 1 regime, is much smoother. For 
smaller values of A, the strong locking behavior of the uncorrelated freely jointed chain characterized through 
a steep slope of the force elongation curve can be captured. 

While the freely jointed chain model of figure|2 left, shows a pronounced locking behavior, the stiffness of the 
wormlike chain model of figureEl right, increases gradually as the locking stretch is approached. For the densely 
packed collagen fibrils considered in the present work, the two parameter wormlike chain model is believed to 
represent the real material behavior more accurately than the single parameter freely jointed chain model which 
might be better suited for randomly oriented polymer chains in rubber. In particular due to the additional 
freedom introduced by the persistence length as a second parameter, we shall thus focus on the wormlike chain 
model for the collagen fibrils in the sequel. 

3 Mechanics of the chain network 

To incorporate the individual chain statistics into an overall constitutive description, we apply an eight chain 
representation of the underlying cooperative macromolecular network structure. These eight chains are em- 
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Figure 3: Kinematics of transversely isotropic eight chain network model 

bedded in a transversely isotropic unit cell with initial cell dimensions a and b. They essentially link at its 
center and extend to the eight individual corners as in figure While the end-to-end length in the undeformed 
configuration is obviously given as ro = i/a 2 + b 2 + b 2 /2, we assume that the end-to-end length in the deformed 

configuration can be expressed as 

r = yjha 2 + [h-h] b 2 /2 (7) 

in terms of the first and fourth invariant I\ and I4. The relevant invariants and their derivatives are given in 
the following form. 



h 
dch 
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h = det(C) 



h = N a 
dch = N a 



(8) 



The first invariant I\ can either be expressed as the trace of the covariant Cauchy Green strain tensor C as 
Ji = G 1 : C or of the contravariant finger tensor b as I\ = g : b. Thereby, the Cauchy Green tensor C is 
defined as the pull back of the covariant spatial metric g whereas the finger tensor b is the push forward of the 
contravariant material metric G _1 . 



C = F ■ g ■ F 



b = F G 



(9) 



Herein, F denotes the deformation gradient as F = Vx¥> with tp being the deformation map between the 
undeformed and the deformed configuration. The determinant of either C or b defines the third invariant -Z3 
which is thus identical to the Jacobian squared I3 — J 2 , whereby J = det (F). The fourth invariant I4 = 
essentially represents the square of the stretch A a along the a direction. As such, I 4 can be expressed as the 
scalar product of the Cauchy Green strain C with the structural tensor N = n ®n , i.e. the dyadic product 
of the normal vectors no in the undeformed configuration. Note that while the stretch A in the out-of-plane 
direction is captured in an affine way through the -Z4 term, the in-plane stretch Af, is obviously not uniquely 
defined. It is thus represented in a non-affine manner through the [I\ — I4] term. The in-plane term [I\ — £4] 
obviously results from the scalar product of C with the remaining contribution G 1 — Nq = G _1 — no <8> no- 
Recall that in the undeformed configuration A a = A& = 1, Ix = 3, I4, = 1 and thus r = \J a 2 + b 2 + b 2 /2. The 
overall energy "3/ of the transversely isotropic eight chain unit cell is assumed to consist of three contributions 



= \]/ blk 4- * chn 



(10) 
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as illustrated in figure|3| The first term ^ blk (/i, J3) captures the effect of bulk incompressibility, e.g. due to a 
surrounding liquid solvent, and is thus of isotropic nature, see Garikipati et al. [19]. The second term $ chn (/i, I4) 
reflects the effective assembly of the individual eight chain energies ip chn . Accordingly, \I/ chn = 7 chn ip chn with 
7 chn denoting the chain density per unit cell. The repulsive term ^ lep (Ii, I4) accounts for a stress- free reference 
configuration and prevents the material from collapsing. The second and third terms essentially depend on 
the key phcnomenological kinematic variable of a single chain, the current end-to-end length r(Ii, I4). For the 
wormlike chain model based on the free energy ip chn = ip wlc according to equation (JHJ, the individual energy 
terms take the following expressions, 
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(11) 
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where we have introduced the following abbreviation for the repulsive weighting factor <3> rep . 

^ = In (if - b2]/2 ) + hn (if ) 



(12) 



The parameter set of the model is thus restricted to the chain density 7 blk , the two wormlike chain parameters, 
i.e. the persistence length A and the contour length L, the cell dimensions a and b and the two bulk parameters 
7 blk and (3. The above introduced free energy '5 defines the Kirchhoff stress r 
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(13) 



as the contravariant push forward of the second Piola Kirchhoff stress 2 dc $ as t = F ■ 2 Ac ^ ■ F t . The 
individual stress contributions follow straightforwardly from the corresponding energy terms introduced in 
equations (fTT|) 
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where g 1 denotes the contravariant spatial metric. In the above equations, we have made use of the following 
abbreviations for the bases of the chain stress and of the repulsive stress. 



x chn = [a 2 -b 2 }N+ b 2 b 

f«p = 1 [a 2 -b 2 ] N + 1 b 2 b. 
h h 



(15) 



Here, N = F ■ No ■ denotes the push forward of the structural tensor TVo. As such, it can be expressed as 
N = n (g> fi in terms of the cell orientation of the deformed configuration fi a — F ■ nP a . In equation (|15l) , we 
have introduced the abbreviation -f chn = F ■ 8 r dcr ■ F^ 1 for the derivative of the end-to-end length r. Since we 



assume a stress- free initial state at r = r Q , we conclude that r rep (ro) = 



_chn 



(r ) and thus r rep (r ) = x chn (r ). 



Accordingly, the repulsive energy contribution 3' rcp introduced in equation (|llf> has been constructed in such 
a way that r rop = F ■ 2dc ^' rop • F t . With the help of the above equations, the spatial Kirchhoff tangent C 
relating the Lie derivative of the Kirchhoff stress LtT to the Lie derivative of the covariant spatial metric Ltg 
as LtT = C : Ltg/2 with 

C = C blk + C chn + C rep (16) 
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is defined through the contravariant push forward of the material tangent 4dc^g)C x I , a $ C 
[F t ®F t ]. Its individual contributions take the following representation. 



[F®F] : 4d c ®c* : 
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The fourth order basis of the chain term C 
can be expressed as follows. 
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In equation i|17|) . 2" denotes the fourth order identity which can be expressed as X = [g^ 1 ® g^ 1 +g~ x ® g 1 ] / 2. 
Here, we have applied the abbreviations Cg> and <g> for the non-standard dyadic products according to the following 
component-wise definitions {•&>o} i j kl = {u}^ <S> {°}ji and {•®o} i:)fe ; = l*}^ {°}jfe- Recall that transverse 
isotropy is basically defined through the normal no of the preferred material direction. This normal will be used 
later on to specify biological remodeling as it is assumed to represent the orientation of collagen fibers within 
a tissue. 

Remark 2 (Special case of an isotropic network model) The classical non-affine isotropic eight 
chain model of Arruda & Boyce [2, 6, 7] which was originally introduced in the context of rubber elasticity 
can be understood as a special case of the present framework. Its undeformed unit cell represents a cube with 
a = b. The chain extension thus reduces to a function of the the first strain invariant 1\ such that r — \JT\aj1. 
For the isotropic eight chain model, the undeformed reference configuration is characterized through A = Ah = 1, 
I\ = 3 and thus ro — V3a/2. The corresponding repulsive term ^ ,rcp = 3/2 In (If ) introduces the stress con- 
tributions T chn — a 2 b and f chn = 3//i a 2 b. Accordingly, the bases for the chain and the repulsive tangent term 
which follow as C chn — a 4 b <g> b and C lop = —6 a 2 /I 2 b <g> b clearly reflect the isotropic nature of this particular 
specification of the model. 



Remark 3 (Special case of a transversely isotropic model) Another special case of the transversely 
isotropic chain network model follows from assuming a degenerated unit cell for which the in-plane dimension 
tends to zero as b = 0. This affine chain model is based on the introduction of chains which are all oriented 
in a single direction Uq. The deformed end-to-end length r = / 2 thus degenerates to a mere function of 

the fourth invariant I4, or rather of the stretch X a , and does no longer incorporate cross-linking network effects. 
In the undeformed case with J4 = 1, the end-to-end length is ro = a/ 2. The repulsive term of equation 



a 4 /2. 



The corresponding bases of the stress and tangent contributions are then 



ichn 



thus degenerates to ^ rep = In (J, 

given exclusively in terms of the structural tensor N as f chn — a 2 N , x rep = 1 / J4 a 2 N , C^"" = a 4 TV <£> TV and 
C rep = —2 //4 a 2 N <S> TV. Obviously, this specific representation of the model, which has been termed "decoupled 
reinforcement model" by Merodio and Ogden (2002), (2003) does not include the characteristic cross-link effects 
of the network structure. 



3.1 Example: Influence of anisotropy 

In the present theory, transverse isotropy is represented by the particular initial cell orientation no and by the 
unit cell dimensions a and b. The influence of both will be studied in the sequel for the illustrative homogeneous 
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load case of uniaxial tension. The eight chains of the model are of wormlike chain type as introduced in section 
12.21 In anticipation of later examples, the contour length L and the persistence length A have been chosen as 
L = 2.125 and A — 1.82. The chain density is taken as 7 chn — 7 x 10 21 . The bulk material of the unit cell is 
defined by the parameters 7 blk = 100 and j3 = 4.5. 



uniaxial tension j10 . uniaxial tension 




stretch stretch 



Figure 4: Influence of anisotropy - variation of cell dimensions and fiber load angle 

Figure left, shows the influence of the cell dimensions for different unit cell heights a with b = 1.95 fixed. 
Thereby the loading axis F is aligned with the fiber direction no- As the cell height increases from a = 1.95, 
i.e. a : b = 1.0, right curve, to a = 2.925, i.e. a : b = 1.5, left curve, the material stiffens considerably. For this 
particular choice of parameters, the material with an aspect ratio of a : b = 1.0, i.e. the isotropic material, has 
a locking stretch which is slightly larger than A* = 1.98 while the locking stress of the material with a : b = 1.5 
is close to A* = 1.24. Recall that at fixed contour length L, fixed persistence length A and fixed cell dimension 
b changes in the cell height a imply changes in the initial end-to-end length r which obviously lead to changes 
in the overall stress strain response. The anisotropic material response is thus highly sensitive to changes in the 
cell dimensions. 

Figure right, illustrates the influence of the orientation of the fiber direction no with respect to the loading 
axis F. For a fixed aspect ratio of a : b — 1.25 and cell dimensions a — 2.43 and b = 1.95, different fiber load 
angles have been studied varying from a = 0°, left curve, to a — 90°, right curve. As expected, the material 
stiffens as the fiber direction rotates towards the loading axis. For a = 90°, the material behaves most compliant 
with a locking stretch of about A* = 1.93, while the locking stretch of the stiffest response at a = 0° is close to 
A* = 1.55. 

3.2 Example: Anisotropic response of rabbit skin 

Finally, the transversely isotropic eight chain model will be applied to simulate the behavior of skin. The 
following simulation is based on an experiment carried out by Lanir & Fung [32], see also Fung [18] for further 
details and Bischoff et al. [5] for corresponding orthotropic eight chain model simulations. In the experiment, 
samples of rabbit skin have been tested parallel and orthogonal to the head to tail direction. Since the collagen 
fibers in skin are basically oriented along the head to tail axis, the skin samples' responses were much stiffer in 
the head to tail direction than orthogonal to it. 

With the set of parameters introduced in the previous section, L — 2.125, A — 1.82, a — 2.43, b — 1.95, 
<y chn — 7 x 10 21 , 7 blk = 100 and j3 = 4.5, the transversely isotropic wormlike chain model nicely captures the 
experimental results by Lanir & Fung [32] which are given in figure left, whereby the data points have been 
reproduced from the original paper. Not only the different locking stretches of about A* = 1.55 and A* = 1.93 
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Figure 5: Anisotropy of rabbit skin - experiment vs. simulation 

are captured nicely by the simulation of figure El right, also the overall shape of the curves corresponds to the 
experimental findings. 



Remark 4 (Freely jointed chain network model) Recall that the use of the freely jointed chain of sec- 
tion 12.11 within the eight chain network model typically overpredicts the locking behavior. Lacking the freedom 
of the second characteristic parameter, the freely jointed chain model experiences difficulties in simultaneously 
adjusting the appropriate locking stretch and the shape of the stress stretch curve. The two-parameter worm- 
like chain model, however, nicely captures both characteristics, not only in the single chain case but also when 
embedded in the representative chain network. 

4 Remodeling of the transversely isotropic tissue 

Soft biological tissues such as muscles, tendons or ligaments show a pronounced orientation of collagen fibers 
along one or two particular directions. Nevertheless, this phenomenon which is typically not observable in 
neonatal tissue only develops upon mechanical loading and is often referred to as functional adaptation. In 
the present section, we suggest a theory that allows for a continuous redistribution of the material's principal 
directions. For the transversely isotropic chain network model, we thus allow the characteristic cell axis to 
rotate according to a particular mechanical stimulus, in our case, the maximum principal strain. 

4.1 Continuum model of remodeling 

The fundamental assumption of the original isotropic eight-chain model by Arruda & Boyce [2,6,7] is that once 
a particular deformation is applied, the unit cell is assumed to rotate instantaneously towards the principal axes. 
The basic idea of the present model is that the unit cell axis no is allowed to gradually align with eigenvector 
n max Q £ c auc hv Green tensor C = F t ■ g ■ F which is associated with maximum eigenvalue A max = max(Ai). 
Here, n™ ax follows straightforwardly from the spectral decomposition of C. 

C = \n\®n\ f = 1,2,3 (19) 

In some biologically relevant cases, we might encounter multiple maximum eigenvalues Aj. For the sake of clarity, 
for the time being, we assume that no alignment takes place for multiple maximum eigenvalues. Following the 
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approach suggested recently by Menzel [35], we introduce the rotation vector u) 

u = — n x n' A nax (20) 

as the scaled vector product of the unit cell orientation no and the direction of the maximum principal strain 
n max ag s k e t c h ec [ m fig ur e |SJ The decomposition of the rotation vector u> into the unit normal n w with 
n u ■ n u = 1 and the magnitude u as 

n x n" lax ||n xn™ ax || , . 

UJ=LUn UJ U W = — * rr uj = — * — - 21 

||n xn^|| t u v ; 

illustrates, that the magnitude of rotation is obviously governed by both the time relaxation parameter t w and 
the angle between no and n™ ax . With the above considerations at hand, the evolution of the unit cell axis no 
can be expressed in the following abstract form 

d t n = u> x n = — [ e • u> ] • n (22) 

whereby e denotes the third order permutation symbol. With the help of equation l|2U|) the above equation can 
be reformulated in the following, maybe more illustrative format. 

dt n = — [ nr x - [ nr* • «o ] "o ] (23) 

Note that for the particular evolution equation (|22() . the orthogonality condition d t no-no = is valid throughout 
the remodeling history. 

4.2 Discrete model of remodeling 

In the present section, we suggest a discrete numerical solution strategy to solve the transient equation (|23|l 
governing the remodeling process. For its temporal discretization, consider a partition of the time interval 
of interest T into a finite number of subintervals T = U/ST 1 [^ fe ,^ fc+1 ] and focus on the typical subinterval 
[t k , t k+1 ] with At = t k+1 — t k denoting the corresponding time increment. In doing so, we assume that the unit 
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cell orientation no is known at t k . The suggested constitutive integrator of the evolution equation (|23(1 is based 
on an exponential integration scheme which allows for a closed form representation of Euler-Rodrigues type. 

rio +1 = exp(- At e • w) • = R ■ (24) 

In the above equation, we have introduced the proper orthogonal tensor R(uj) = R(u>, n u ) which is characterized 
through the following closed form expression. 

R = cos( At a) I - sin( At u) e ■ n u + [ 1 - cos( At w) ] n u <g n w (25) 

The combination of equations l|24|l and (|25|l straightforwardly renders the following update formula for the cell 
axis n T , 

riQ +1 = cos( At w)nj + sin (At w) x n\ + [ 1 - cos(At us) ] [ n" • ] n B (26) 

as in Menzel [35]. Recall from equation H20|) . that the current rotation vector uj = uj (n|j +1 ) is a function of the 
new direction ng +1 . Accordingly, the above equation is actually an implicit update equation in the unit cell 
direction ng +1 with dependencies on rig +1 on the left and righthand side. In what follows, however, we shall 
tacitly assume that n.g +1 x n w w nj x n", such that the rotation vector u) can be approximated as a function 
of the known orientation rig . Accordingly, u) can then be updated explicitly as 

uj = — n£ x n^ ax (27) 

which is a reasonable assumption in the context of the gradual alignment postulated herein. 



Remark 5 (Classical eight chain model) Recall that the unit cell edges of the classical isotropic eight 
chain model of Arruda & Boyce [2,6,7] are assumed to align instantaneously with the axes of principal strain. For 
isotropic rubber elasticity, this might be a reasonable assumption. However, it is widely accepted that biological 
tissues show a gradual alignment of the material's principal axes with respect to a mechanical stimulus due to 
a cascade of cellular and molecular events involved in the synthesis and breakdown of collagen that occurs over 
time. The suggested remodeling approach inherently accounts for the successive reorientation of the anisotropic 
unit cell with respect to the direction of maximum principal strain. Accordingly, the final biological equilibrium 
state is always characterized through the alignment of the unit cell axes rig with the maximum principal strain 



4.3 Example: Influence of remodeling 

The features of the suggested remodeling approach are demonstrated for the homogeneous problem under 
uniaxial tension illustrated in figure The model parameters are adopted from the previous examples of 
section|21 In addition, the relaxation parameter is chosen to r = 1.0 and the time step size is At = 0.1. Similar 
to the previous examples, the specimen is loaded up to a final load of F — 20000 in the first ten load steps. 
However, now, the load is held constant for another 90 time steps to allow for a smooth adaptation of the cell 
orientation. 

Starting at a fiber load angle of a = 90°, the cell axis no gradually rotates towards the loading axis, i.e. a = 0°, 
as time evolves, as illustrated in figure left. The curves nicely demonstrate the magnitude of adaptation 
uj = | |ng x n™ ax | \/t u which is obviously large for a large mismatch of axes and which decreases as the fiber load 
angle tends to zero. The additional influence of the relaxation parameter t u manifests itself in a fast adaptation 
for small values of t u , e.g. for r w = 0.5, left curve. For larger values of t u , e.g. for t w = 8.0, right curve, the 
adaptation time increases. The suggested model is thus able to capture a gradual reorientation of the axis of 
transverse isotropy with respect to the principal strain axis. 

The local strengthening upon reorientation is nicely documented by figured right. It shows the evolution of 
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the stretch as a function of time. In the loading phase, the stretch increases up to about A* = 1.93. This value 
corresponds to the locking stretch of the a = 90° orientation in figure 0] right. Then, upon remodeling, the 
specimen contracts as the strong material axis rotates towards the loading axis. The plateau indicates the state 
of biological equilibrium. The alignment of the cell axis with the axis of loading reduces the overall stretch to 
A* = 1.55 which agrees nicely with the a = 0° results of figure H right. 

4.4 Example: Fiber reorientation in tendons 

Let us now turn to the simulation of a real biomechanical boundary value problem motivated by an experiment of 
engineered functional tendon constructs as documented by Calve et al. [9] . Experimental observations confirm, 
that in the absence of loading, the in vitro grown tendon constructs showed the typical characteristics of 
embryonic tendon demonstrated by the absence of a collageneous scaffold. It is only upon mechanical loading 
that collagen fibrils form within the tendon and orient themselves along the loading direction. The long term goal 
of the present project is thus to experimentally analyze and computationally predict mechanically stimulated 
remodeling in the form of fiber reorientation. In the future, other microstructural changes such as the increase 
in cross linkage or the post-depositional fusion of fibrils will be considered as well. These might be incorporated 
straightforwardly in the present framework through changes in the contour length or in the persistence length 
for which appropriate evolution equations have to be defined. 

As a preliminary study, we carry out a finite element simulation of the remodeling process in a cylindrical 
model tendon. The tendon material is described with the transversely isotropic wormlike chain model whereby 
an initially random fiber orientation is assumed to represent the neonatal state. The current fiber orientation no 
is introduced as an internal variable which is stored locally on the integration point level. The tendon structure 
has an initial cross section of unit area and a length of twelve units, respectively. It is discretized with 2304 
eight-noded brick elements introducing about 8000 degrees of freedom. The two wormlike chain parameters, i.e. 
the the contour length and the persistence length, take values of L = 2.50 and A = 1.82. Note that all lengths 
in the model are normalized by the segment length I. The bulk parameters take the values of 7 blk = 100 and 
j3 = 4.5, the chain density is chosen to 7 chn = 7 x 10 21 , the absolute temperature is 6 = 310. The dimensions 
of the transversely isotropic unit cell are taken as a = 2.43 and b = 1.85, again normalized with respect to the 
segment length I. The relaxation time for the remodeling procedure is chosen as r=1.0. 

In the first 40 time steps of At — 0.0025, the model tendon is gradually stretched to about 200 % of its initial 
length. The final load of F — 1000 is then held constant for another 40 time steps of At = 0.01 to allow for 
fiber reorientation, see figure |H1 left. Figure |S1 right, shows the temporal evolution of the elongation of the 
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Figure 8: Remodeling of soft biological tissue - prescribed loading and extension 

tendon. As the tendon is loaded, i.e. for 0.0 < t < 1.0, the elongation increases smoothly up to almost 12 units. 
During the reorientation period, i.e. for 1.0 < t < 5.0, the tendon obviously stiffens considerably due to fiber 
reorientation. Accordingly, the elongation reduces to about 8.5 units at the final biological equilibrium state, 
see figure |H1 right. 

Figure [5] depicts six representative stages of the remodeling history. The top figure shows the initial state at the 
beginning of the loading history. The arrows indicate the the initial orientation of the unit cells no, or, in the 
biomechanical sense, the directions of pronounced collagen fiber orientation, that have been assigned randomly 
to each individual integration point. The contour levels ranging from white to black depict the fiber load angle 
a varying from a = 90°, i.e. the strong cell axis being orthogonal to the loading axis, to a = 0°, i.e. a full 
alignment of the cell axis with the axis of loading. 

Figure|nifrom top to bottom nicely documents the history of remodeling. While the tendon is loaded uniaxially, 
a clear reorientation of the collagen fibers with respect to the loading axis can be observed. In this sense, the 
last figure of the series represents a biological equilibrium state, for which dt no = 0. No further remodeling 
takes place as all fibers are aligned with the maximum principal strain direction. Accordingly, the fiber load 
angle indicated through the underlying contour plots evolves gradually from a random distribution in the top 
figure to a uniform distribution with a fiber load angle of a = 0° in the bottom figure. 

While converging towards the final biological equilibrium state the suggested procedure showed a remarkably 
stable algorithmic behavior. Recall that during the loading process, the tendon has been stretched to twice its 
original length. The stretch of the tendon as well as the remarkable reduction of its cross section are clearly 
visible in figure The algorithm is thus able to capture both, kinematic and constitutive non-linearities. 
Apparently, the explicit update of the fiber direction does not lead to computational instabilities as long as the 
relaxation parameter is sufficiently large and the time step size is chosen sufficiently small. 

5 Discussion 

A new transversely isotropic chain network model has been proposed which is particularly suited to simulate 
remodeling of collagen fibers in soft biological tissue. For the individual chains, a wormlike chain model was 
adopted. In contrast to the classical uncorrelated freely jointed chain, the correlated wormlike chain shows an 
initial stiffness such that its curvature changes smoothly in the contour space. The stiffness, or rather the degree 
of correlation, is governed by a second parameter besides the contour length, the so-called persistence length. 
With the additional freedom of this second parameter, the wormlike chain model can capture the behavior of 
the freely jointed chain and a beam-like behavior as limit cases. 
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The chain network is characterized by a representative eight chain unit cell. In contrast to the cubic cell of 
the isotropic eight-chain model, the present unit cell has one preferred direction which characterizes the axis of 
transverse isotropy. As such, the present model can be understood as a generalization of the classical eight chain 
model which captures the original isotropic eight chain model as a special limit case with equal cell dimensions. 
In the other limit, i.e. when the in-plane dimensions of the unit cell tend to zero, the model captures the effects 
of unidirectional fiber reinforcement neglecting cross-link effects of the network structure. 

To incorporate biomechanical remodeling, the characteristic axis of the transversely isotropic unit cell was 
allowed to rotate with respect to a biological stimulus, in our case the maximum principal strain axis. A 
theoretical framework of remodeling and its numerical realization based on an exponential update scheme of 
Euler-Rodrigues type have been introduced. Thereby, the characteristic unit cell axis has been introduced as 
an internal variable on the integration point level of a finite element realization. 

A first comparison of the suggested approach with the experimental findings of in vitro engineered tendon 
constructs revealed a remarkably good agreement. Being essentially based on micromechanical considerations, 
the present model is governed by a limited number of physically- motivated material parameters. As such, it is 
believed to be ideally suited to simulate not only the passive behavior of soft biological tissues but also their 
active response to changes in the mechanical loading environment. 
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